module transform_module

contains

  subroutine coordinate_transform()
    use my_mpi
    use logs
    use camera_module
#ifdef VELOCITY
    use particles    , only : NumPart_ByPE, readCoordinates, readVelocity
#else
    use particles    , only : NumPart_ByPE, readCoordinates
#endif
    implicit none
    integer*4 :: ipart
    real*4                 :: theta_rad,phi_rad,psi_rad
    real*4, dimension(3)   :: pos,tmp
#ifdef VELOCITY
    real*4, dimension(3)   :: vel
#endif
    real*4, dimension(3,3) :: RotMatrix
    real*8, parameter      :: pi = 3.14159265358979323846d0

#ifdef USE_EIGENVECTORS
    
    write(logunit,*)'  -transforming with eigenvectors'
    call flush(logunit)

    ! Just dot product each position with the new basis vectors
    do ipart=1,NumPart_ByPE(MyPE)
       pos(1) = readCoordinates(1,ipart) - x0
       pos(2) = readCoordinates(2,ipart) - y0
       pos(3) = readCoordinates(3,ipart) - z0
       tmp(1) = pos(1)*xprime(1) + pos(2)*xprime(2) + pos(3)*xprime(3)
       tmp(2) = pos(1)*yprime(1) + pos(2)*yprime(2) + pos(3)*yprime(3)
       tmp(3) = pos(1)*zprime(1) + pos(2)*zprime(2) + pos(3)*zprime(3)
       readCoordinates(1,ipart) = tmp(1) + x0
       readCoordinates(2,ipart) = tmp(2) + y0
       readCoordinates(3,ipart) = tmp(3) + z0
#ifdef VELOCITY
       vel(1) = readVelocity(1,ipart) 
       vel(2) = readVelocity(2,ipart) 
       vel(3) = readVelocity(3,ipart) 
       tmp(1) = vel(1)*xprime(1) + vel(2)*xprime(2) + vel(3)*xprime(3)
       tmp(2) = vel(1)*yprime(1) + vel(2)*yprime(2) + vel(3)*yprime(3)
       tmp(3) = vel(1)*zprime(1) + vel(2)*zprime(2) + vel(3)*zprime(3)
       readVelocity(1,ipart) = tmp(1) 
       readVelocity(2,ipart) = tmp(2) 
       readVelocity(3,ipart) = tmp(3) 
#endif       
    enddo
#else

    ! Convert angles to radians
    theta_rad = pi*theta / 1.8e2    
    phi_rad   = pi*phi   / 1.8e2    
    psi_rad   = pi*psi   / 1.8e2    

    ! Construct rotation matrices
    RotMatrix(1,1) =  cos(phi_rad)*cos(psi_rad)
    RotMatrix(1,2) = -cos(phi_rad)*sin(psi_rad) 
    RotMatrix(1,3) =  sin(phi_rad) 
    RotMatrix(2,1) =  cos(theta_rad)*sin(psi_rad) + sin(theta_rad)*sin(phi_rad)*cos(psi_rad)
    RotMatrix(2,2) =  cos(theta_rad)*cos(psi_rad) - sin(theta_rad)*sin(phi_rad)*sin(psi_rad)
    RotMatrix(2,3) = -sin(theta_rad)*cos(phi_rad)
    RotMatrix(3,1) =  sin(theta_rad)*sin(psi_rad) - cos(theta_rad)*sin(phi_rad)*cos(psi_rad)
    RotMatrix(3,2) =  sin(theta_rad)*cos(psi_rad) - cos(theta_rad)*sin(phi_rad)*sin(psi_rad)
    RotMatrix(3,3) =  cos(theta_rad)*cos(phi_rad)

    ! Rotate about each axis with a matrix operation
    do ipart=1,NumPart_ByPE(MyPE)
       readCoordinates(:,ipart) = matmul(RotMatrix,readCoordinates(:,ipart))
#ifdef VELOCITY
       readVelocity(:,ipart)    = matmul(RotMatrix,readVelocity(:,ipart))
#endif
    enddo
#endif

  end subroutine coordinate_transform

end module transform_module




